#### Meta Data ####

# article: "The Ambivalent Effect of Autocratization on Domestic Terrorism" 
# journal: "Studies in Conflict & Terrorism"
# authors: "Lott, Lars; Croissant, Aurel; Trinn. Christoph"
# date: 2023-10-02
# written under "R version 4.3.1 (2023-06-16 ucrt)"

#### Preliminaries ####

R.version$version.string

# clear workspace
rm(list=ls())

# set working directory
getwd()

#### Preliminaries ####

library(MASS)
library(multiwayvcov)
library(ggeffects)
library(PanelMatch)
library(tidyverse)
library(ggpubr)
library(here)
library(haven)
library(imputeTS)
library(modelsummary)



#### Load Dataset ####

vdem <- readRDS("data/vdem12.rds") 

summary(vdem$aut_ep_outcome_agg)

#### Define Treatment ####

vdem <- vdem %>%
  mutate(autocratization_treatment = aut_ep,
         democratic_breakdown = ifelse(aut_ep_outcome_agg == 1, 1, 0), 
         democratic_recession = ifelse(aut_ep_outcome_agg == 2, 1, 0), 
         regressed_autocracy = ifelse(aut_ep_outcome_agg == 3, 1, 0))

table(vdem$autocratization_treatment)
table(vdem$democratic_breakdown)
table(vdem$democratic_recession)
table(vdem$regressed_autocracy)

vdem <- vdem %>%
  mutate(democratic_breakdown_recession = ifelse(democratic_breakdown == 1 |democratic_recession == 1, 1,0 ))

table(vdem$democratic_breakdown_recession)



vdem <- vdem %>%
  mutate(log_count_attacks_dom = log(count_attacks_dom +0.01),
         log_count_attacks_int = log(count_attacks_int +0.01), 
         log_count_attacks = log(count_attacks +0.01),
         log_sum_nkill_dom = log(sum_nkill_dom +0.01),
         log_sum_nkill_int = log(sum_nkill_int +0.01),
         log_sum_nkill = log(sum_nkill +0.01)) %>%
  filter(year >= 1966)


vdem <- vdem %>%
  group_by(country_id) %>%
  mutate(log_count_attacks_dom = na_interpolation(log_count_attacks_dom, option = "spline",  maxgap = 1), 
         log_sum_nkill_dom = na_interpolation(log_sum_nkill_dom, option = "spline",  maxgap = 1), 
         count_attacks_dom = na_interpolation(count_attacks_dom, option = "spline",  maxgap = 1), 
         sum_nkill_dom = na_interpolation(sum_nkill_dom, option = "spline",  maxgap = 1), 
         log_count_attacks_dom = ifelse(log_count_attacks_dom<0, 0, log_count_attacks_dom),
         log_sum_nkill_dom = ifelse(log_sum_nkill_dom<0, 0, log_sum_nkill_dom),
         count_attacks_dom = ifelse(count_attacks_dom<0, 0, count_attacks_dom),
         sum_nkill_dom = ifelse(sum_nkill_dom<0, 0, sum_nkill_dom))

vdem <- vdem %>%
  dplyr::select(country_id, country_name, year, aut_ep, 
                aut_ep_id,  log_count_attacks_dom, log_count_attacks_int, log_count_attacks, 
                log_sum_nkill_dom, log_sum_nkill_int, log_sum_nkill, count_attacks_dom, count_attacks_int, count_attacks, 
                sum_nkill_dom, sum_nkill_int, sum_nkill,e_gdppc, e_pop,
                v2x_polyarchy, milper_extrapol, eprratio_extrapol, democratic_recession, democratic_breakdown, 
                civil_war_ucdp, v2x_regime, aut_ep_outcome, aut_ep_outcome_agg, 
                v2caviol, dem_ep, dem_ep_id, dem_ep_outcome_agg, dem_ep_outcome, e_regionpol, e_regionpol_6C, 
                autocratization_treatment, regressed_autocracy, democratic_breakdown_recession, ma_count_attacks_dom, ma_sum_nkill_dom, democracy_stock, 
                v2x_veracc, v2x_diagacc, v2x_horacc, e_pelifeex, e_peaveduc, e_area) %>%
  filter(year <=2020) %>%
  group_by(country_id) %>%
  fill(e_area)

vdem <- vdem %>%
  drop_na(autocratization_treatment) %>%
  mutate(e_gdppc_sqaured = e_gdppc^2, 
         e_area_log = log10(e_area))

summary(vdem$e_area_log)
summary(vdem$e_gdppc_sqaured)

#------------------------------------------------------------------------------#
#------------------------------ Figure D1 -------------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_10 <- c(0:10)


vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$democratic_breakdown_recession)
table(vdem$democratic_breakdown_recession)


#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "democratic_breakdown_recession", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_10, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "democratic_breakdown_recession", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_10, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "democratic_breakdown_recession", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_10, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) +  I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) +  
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_10, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 25 democratic breakdown

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "democratic_breakdown_recession", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_D1/Figure_Frequency_Distribution_10_years.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for democratic breakdown/ recession", rot = 90))
ggsave("results/Figure_D1/Figure_Covariate_Balance_Autocratization_10_years.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_D1/Figure_Improved_Covariate_Balance_10_years.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "democratic_breakdown_recession", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####
 
#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5,6,7,8,9,10))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5,6,7,8,9,10))


saveRDS(M1_pred, "results/Figure_D1/data_M1_pred_dem_breakdown_10_years.rds")
saveRDS(M2_pred, "results/Figure_D1/data_M2_pred_dem_breakdown_10_years.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5,6,7,8,9,10))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5,6,7,8,9,10))

saveRDS(M1_pred_ci09, "results/Figure_D1/data_M1_pred_dem_breakdown_ci09_10_years.rds")
saveRDS(M2_pred_ci09, "results/Figure_D1/data_M2_pred_dem_breakdown_ci09_10_years.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_D1/data_M1_pred_dem_breakdown_10_years.rds") 
data_figure_1b <- readRDS("results/Figure_D1/data_M2_pred_dem_breakdown_10_years.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_D1/data_M1_pred_dem_breakdown_ci09_10_years.rds") 
data_figure_1b_09 <- readRDS("results/Figure_D1/data_M2_pred_dem_breakdown_ci09_10_years.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Democratic Breakdown/Recession", 
       y = "Estimated effect of Democratic Breakdown/Recession", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a 

ggsave("outputs/Figure_D1.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_D1.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#---------------------------- Figure D2   -------------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_10 <- c(0:10)

vdem$country_id <- as.integer(vdem$country_id)

summary(vdem$regressed_autocracy)
table(vdem$regressed_autocracy)

#### 1. Effect of Autocratization on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "regressed_autocracy", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_10, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "regressed_autocracy", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_10, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "regressed_autocracy", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_10, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_10, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 46 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "regressed_autocracy", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_D2/Figure_Frequency_Distribution_10_years.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for Autocratization", rot = 90))
ggsave("results/Figure_D2/Figure_Covariate_Balance_Autocratization_10_years.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom", "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_D2/Figure_Improved_Covariate_Balance_10_years.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "regressed_autocracy", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5,6,7,8,9,10))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5,6,7,8,9,10))


saveRDS(M1_pred, "results/Figure_D2/data_M1_pred_dem_breakdown_10_years.rds")
saveRDS(M2_pred, "results/Figure_D2/data_M2_pred_dem_breakdown_10_years.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5,6,7,8,9,10))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time =c(0,1,2,3,4,5,6,7,8,9,10))

saveRDS(M1_pred_ci09, "results/Figure_D2/data_M1_pred_dem_breakdown_ci09_10_years.rds")
saveRDS(M2_pred_ci09, "results/Figure_D2/data_M2_pred_dem_breakdown_ci09_10_years.rds")


#### Regressed Autocracy Effect ####

data_figure_1a <- readRDS("results/Figure_D2/data_M1_pred_dem_breakdown_10_years.rds") 
data_figure_1b <- readRDS("results/Figure_D2/data_M2_pred_dem_breakdown_10_years.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_D2/data_M1_pred_dem_breakdown_ci09_10_years.rds") 
data_figure_1b_09 <- readRDS("results/Figure_D2/data_M2_pred_dem_breakdown_ci09_10_years.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocratic Consolidation", 
       y = "Estimated effect of Autocratic Consolidation", 
       title = "Propensity Score Weighting: Terrorist Attacks")
f1a

ggsave("outputs/Figure_D2.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_D2.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)


#------------------------------------------------------------------------------#
#------------------------------------ Figure D3 -------------------------------#
#------------------------------------------------------------------------------#

vdem <- distinct(vdem, country_id, year, .keep_all= TRUE) # Distinct observations

vdem$year <- as.integer(vdem$year)


vdem <- as.data.frame(vdem)
class(vdem)

# Define Lags and leads for all

lag_2 <- c(1:2)
lead_10 <- c(0:10)

vdem <- vdem %>%
  dplyr::select(-c(country_name, aut_ep_id))

vdem$country_id <- as.integer(vdem$country_id)

#### 1. Effect of democratic_breakdown_delta on Terrorism in all regimes #####

## Before Refinement ##
PM.results.none <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                              treatment = "autocratization_treatment", refinement.method = "none", 
                              data = vdem, match.missing = FALSE, 
                              covs.formula = ~ I(lag(ma_count_attacks_dom, lag_2)), 
                              size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                              lead = lead_10, forbid.treatment.reversal = FALSE, 
                              use.diagonal.variance.matrix = TRUE)

## Mahanlanobis Distance Matching ##
PM.results_mahalanobis <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                     treatment = "autocratization_treatment", refinement.method = "mahalanobis", 
                                     data = vdem, match.missing = FALSE, 
                                     covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                       I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                       I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                       I(lag(ma_count_attacks_dom, lag_2)), 
                                     size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                     lead = lead_10, forbid.treatment.reversal = FALSE, 
                                     use.diagonal.variance.matrix = TRUE)

## Propensity Score Matching ##
PM.results_psmatch <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                 treatment = "autocratization_treatment", refinement.method = "ps.match", 
                                 data = vdem, match.missing = FALSE, 
                                 covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                   I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                   I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                   I(lag(ma_count_attacks_dom, lag_2)), 
                                 size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                 lead = lead_10, forbid.treatment.reversal = FALSE,
                                 use.diagonal.variance.matrix = TRUE)

## Propensity Score Weighting ##
PM.results_psweight <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                                  treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                                  data = vdem, match.missing = FALSE, 
                                  covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                                    I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                                    I(lag(e_pop, lag_2)) + I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                                    I(lag(ma_count_attacks_dom, lag_2)), 
                                  size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                                  lead = lead_10, forbid.treatment.reversal = FALSE,
                                  use.diagonal.variance.matrix = TRUE)

class(PM.results_psweight)

msets <- PM.results_psweight$att

print(msets) # 71 treatments

######################################################

## Display different matched sets ##
mset <- PM.results_psweight$att[10]

DisplayTreatment(unit.id = "country_id",
                 time.id = "year", legend.position = "none",
                 xlab = "year", ylab = "Country Code",
                 treatment = "autocratization_treatment", data = vdem,
                 matched.set = mset, # this way we highlight the particular set
                 show.set.only = TRUE)

summary(PM.results_mahalanobis$att)
summary(PM.results_psweight$att)


par(mfrow=c(1,1)) 
plot(PM.results_psweight$att, main ="Distribution of Matched Set Sizes")

dev.copy(pdf,'results/Figure_D3/Figure_Frequency_Distribution_10_years.pdf', width = 10, height = 5)
dev.off()


#### Check Improved Covariate Balance due to Matching over the Pre-Treatment Time ####

## Get covariate balance ##

df_none <- get_covariate_balance(PM.results.none$att,
                                 data = vdem,
                                 covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                 plot = FALSE)

df_mahalanobis <- get_covariate_balance(PM.results_mahalanobis$att,
                                        data = vdem,
                                        covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                       "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                        plot = FALSE)


df_psmatch <- get_covariate_balance(PM.results_psmatch$att,
                                    data = vdem,
                                    covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                   "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                    plot = FALSE)


df_psweight <- get_covariate_balance(PM.results_psweight$att,
                                     data = vdem,
                                     covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                                                    "e_pop", "ma_count_attacks_dom", "democracy_stock", "e_area_log"),
                                     plot = FALSE)


#Tidying outputs #
df_none <- as.data.frame(df_none)

df_none <- df_none %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F1 <- ggplot(df_none, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Before Refinement") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_mahalanobis <- as.data.frame(df_mahalanobis)

df_mahalanobis <- df_mahalanobis %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F2 <- ggplot(df_mahalanobis, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Mahanlanobis Distance\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

df_psmatch <- as.data.frame(df_psmatch)

df_psmatch <- df_psmatch %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F3 <- ggplot(df_psmatch, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Matching") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()


df_psweight <- as.data.frame(df_psweight)

df_psweight <- df_psweight %>%
  rownames_to_column("year") %>%
  pivot_longer(-year, names_to = "covariate", values_to = "sd_covariate") %>%
  mutate(year = case_when(year =="t_4" ~ -4, 
                          year =="t_3" ~ -3,
                          year =="t_2" ~ -2,
                          year =="t_1" ~ -1,
                          year == "t_0" ~ 0))

F4 <- ggplot(df_psweight, aes(x = year, y = sd_covariate, color = covariate)) +
  geom_line() +
  geom_hline(yintercept=0, color = "black") +
  labs(x = "", y = "", 
       title = "Propensity Score\n Weighting") +
  ylim(-0.8, 0.8) +
  theme(legend.position = "none", 
        plot.title = element_text(hjust = 0.5)) +
  theme_pubr()

Figure1 <- ggarrange(F1, F2, F3, F4,
                     ncol = 4, nrow = 1, 
                     common.legend = TRUE)

annotate_figure(Figure1, 
                left = text_grob("Standardized Mean Differences\n for autocratization", rot = 90))
ggsave("results/Figure_D3/Figure_Covariate_Balance_Autocratization_10_years.pdf", height = 15, width = 35, units= c("cm"))


#### Improved Covariate Balance due to the Refinement of Matched Sets. 

par(mfrow=c(1,3)) 
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_mahalanobis$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psmatch$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))
balance_scatter(non_refined_set = PM.results.none$att,
                refined_list = list(PM.results_psweight$att),
                data = vdem,
                covariates = c("e_gdppc",  "milper_extrapol", "eprratio_extrapol",
                               "e_pop",  "ma_count_attacks_dom",  "democracy_stock", "e_area_log"))

dev.copy(pdf,'results/Figure_D3/Figure_Improved_Covariate_Balance_10_years.pdf', width = 10, height = 5)
dev.off()

##### Democratic Breakdown as Treatment and Terrorist Attacks ####

## domestic terrorist attacks ##

PM_M1_aut <- PanelMatch(lag = 2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + 
                          I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_count_attacks_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "count_attacks_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)

M1_aut<- PanelEstimate(sets = PM_M1_aut, data = vdem)

M1_aut_ci09 <- PanelEstimate(sets = PM_M1_aut, data = vdem, confidence.level = .9)


## terrorist fatalities domestic ##

PM_M2_aut <- PanelMatch(lag =2, time.id = "year", unit.id = "country_id", 
                        treatment = "autocratization_treatment", refinement.method = "ps.weight", 
                        data = vdem, match.missing = FALSE, 
                        covs.formula = ~ I(lag(e_gdppc, lag_2)) + 
                          I(lag(milper_extrapol, lag_2)) + I(lag(eprratio_extrapol, lag_2)) + 
                          I(lag(e_pop, lag_2)) +   I(lag(democracy_stock, lag_2)) + I(lag(e_area_log, lag_2)) + 
                          I(lag(ma_sum_nkill_dom, lag_2)), 
                        size.match = 10, qoi = "att" ,outcome.var = "sum_nkill_dom",
                        lead = lead_10, forbid.treatment.reversal = FALSE,
                        use.diagonal.variance.matrix = TRUE)
M2_aut <- PanelEstimate(sets = PM_M2_aut, data = vdem)
M2_aut_ci09 <- PanelEstimate(sets = PM_M2_aut, data = vdem, confidence.level = .9)

par(mfrow=c(1,2)) 
plot(M1_aut)
plot(M2_aut)

#### Save results as data.frames to build coherent predicted plots with 2WFE and PanelMatch Models ####

#95% CIs ##

M1_pred <- summary(M1_aut)

M1_pred <- M1_pred$summary

M1_pred <- cbind(M1_pred , time = c(0,1,2,3,4,5, 6, 7, 8, 9, 10))

M2_pred <- summary(M2_aut)

M2_pred <- M2_pred$summary

M2_pred <- cbind(M2_pred , time = c(0,1,2,3,4,5,6, 7, 8, 9, 10))


saveRDS(M1_pred, "results/Figure_D3/data_M1_pred_dem_breakdown_10_years.rds")
saveRDS(M2_pred, "results/Figure_D3/data_M2_pred_dem_breakdown_10_years.rds")

# 90% CIs ##

M1_pred_ci09 <- summary(M1_aut_ci09)

M1_pred_ci09 <- M1_pred_ci09$summary

M1_pred_ci09 <- cbind(M1_pred_ci09 , time = c(0,1,2,3,4,5,6, 7, 8, 9, 10))

M2_pred_ci09 <- summary(M2_aut_ci09)

M2_pred_ci09 <- M2_pred_ci09$summary

M2_pred_ci09 <- cbind(M2_pred_ci09 , time = c(0,1,2,3,4,5,6, 7, 8, 9, 10))

saveRDS(M1_pred_ci09, "results/Figure_D3/data_M1_pred_dem_breakdown_ci09_10_years.rds")
saveRDS(M2_pred_ci09, "results/Figure_D3/data_M2_pred_dem_breakdown_ci09_10_years.rds")


#### Autocratization Effect ####

data_figure_1a <- readRDS("results/Figure_D3/data_M1_pred_dem_breakdown_10_years.rds") 
data_figure_1b <- readRDS("results/Figure_D3/data_M2_pred_dem_breakdown_10_years.rds") 

data_figure_1a <- as.data.frame(data_figure_1a)
data_figure_1b <- as.data.frame(data_figure_1b)

data_figure_1a_09 <- readRDS("results/Figure_D3/data_M1_pred_dem_breakdown_ci09_10_years.rds") 
data_figure_1b_09 <- readRDS("results/Figure_D3/data_M2_pred_dem_breakdown_ci09_10_years.rds") 

data_figure_1a_09 <- as.data.frame(data_figure_1a_09)
data_figure_1b_09 <- as.data.frame(data_figure_1b_09)

data_figure_1a <- cbind(data_figure_1a, data_figure_1a_09[,3:4])
data_figure_1b <- cbind(data_figure_1b, data_figure_1b_09[,3:4])


#### Building Subplots for each Estimand ####

f1a <- ggplot(data = data_figure_1a, aes(x = time)) +
  geom_linerange(aes(ymin=`2.5%`, ymax = `97.5%`), size = 1.3, alpha = 0.7, color = "#a6a6a6") +
  geom_linerange(aes(ymin=`5%`, ymax = `95%`), size = 1.5, color = "#5a5a5a") +
  geom_point(aes(y=estimate), size = 4, color = "black") +
  geom_hline(aes(yintercept = 0), color = "red", linetype = "dashed", alpha = 0.5) +
  theme_pubr()+
  scale_x_continuous(breaks=seq(0,10,1)) +
  labs(x = "Year around Autocraitzation", 
       y = "Estimated effect of Autocraitzation", 
       title = "Propensity Score Weighting: Terrorist Attacks")

f1a
ggsave("outputs/Figure_D3.png", units= c("cm"), width = 20, height = 12, dpi = 900)
ggsave("outputs/Figure_D3.pdf", units= c("cm"), width = 20, height = 12, dpi = 900)



